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Non-Oscillatory  Central  Schemes 
for  One-  and  Two-Dimensional  MHD  Equations.  I.* 

Jorge  Baibas}  Eitan  Tadmoh  and  Cheng-Chin  Win 
August  20,  2003 


Abstract 

In  this  paper  we  utilize  a  family  of  high-resolution,  non-oscillatory  central  schemes  for  the  approximate 
solution  of  the  equations  of  ideal  magnetohydrodynamics  (MHD)  in  one-  and  two-space  dimensions.  We 
present  several  prototype  problems.  Solutions  of  one-dimensional  shock-tube  problems  is  carried  out  using 
second-  and  third-order  central  schemes  [19,  18],  and  we  use  the  second-order  central  scheme  [11]  which  is 
adapted  for  the  solution  of  the  two-dimensional  Kelvin-Helmholtz  and  Orszag-Tang  problems.  A  qualitative 
comparison  reveals  an  excellent  agreement  with  previous  results  based  on  upwind  schemes.  Central  schemes, 
however,  require  little  knowledge  about  the  eigen-structure  of  the  problem  —  in  fact,  we  even  avoid  an  explicit 
evaluation  of  the  corresponding  Jacobians,  while  at  the  same  time  they  eliminate  the  need  for  dimensional 
splitting.  The  one-  and  two-dimensional  computations  reported  in  this  paper  demonstrate  the  remarkable 
versatility  of  central  schemes  as  black-box,  Jacobian-free  MHD  solvers. 

AMS  subject  classification:  Primary  65M10;  Secondary  65M05 

Key  words.  Multidimensional  conservation  laws,  ideal  Magnetohydrodynamics  (MHD)  equations,  high-resolution 
central  schemes,  non-oscillatory  reconstructions,  Jacobian-free  form. 


1  Introduction 

In  this  paper  we  present  second-  and  third-order  non-oscillatory  central  schemes  for  the  approximate  solution 
of  the  equations  of  ideal  Magnetohydrodynamics 


Pt 

=  -V  •  (pv), 

(1.1) 

0v)t 

=  —  V  •  [pvvT  +  (p  +  i-B2)/3X3  —  BBt], 

(1.2) 

=  V  x  (• 

v  x  B), 

(1.3) 

et 

=  -V  • 

/-y  If) 

( - p  -\ — pv2)v  —  (v  x  B)  x  B 

Lv7  —  1  2  '  J 

(1.4) 

Here,  p  and  e  are  scalar  quantities  representing  respectively,  the  mass  density  and  the  total  internal  energy, 
v  =  [vx,vy,vz)T  is  the  velocity  field  with  L2-norm  v2  :=  ||v||2,  and  B  =  (Bx,  Byi  BZ)T  and  B2  :=  ||B||2 
represent  the  magnetic  field  and  its  L-norm.  Finally,  the  pressure  p  is  coupled  to  the  internal  energy,  e  = 
\pv2  +  ^B2  +p/ (7—  1),  where  7  is  the  (fixed)  ratio  of  specific  heats.  The  system  is  augmented  by  the  solenoidal 
constraint  V  •  B  =  0;  that  is,  if  the  condition  V  •  B  =  0  is  satisfied  initialy  at  t  =  0,  then  by  (1.3)  it  remains 
invariant  in  time. 

The  intrinsic  complexity  of  these  MHD  equations  suggests  the  class  of  central  schemes  as  an  efficient  alter¬ 
native  for  the  class  of  upwind  schemes,  for  computing  approximate  solutions  the  MHD  problems  (1.1)-(1.4). 

‘Research  was  supported  in  part  by  NSF  Grant  No.  01-07428  (J.B.  and  E.T.),  by  ONR  Grant  No.  N00014-91-J-1076  (E.T.) 
NSF  VIGRE  Grant  (J.B.)  and  by  NASA  grant  No.  NAG5-12986  (C.C.  W.). 
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t  Department  of  Mathematics,  Center  of  Scientific  Computation  And  Mathematical  Modeling  (CSCAMM)  and  Institute  for 
Physical  Science  and  Technology  (IPST),  University  of  Maryland,  MD  20742.  e-mail:  tadmor@cscamm.umd.edu 

§  Department  of  Physics  and  Astronomy,  University  of  California,  Los  Angeles,  CA  90095,  email:  wu@physics.ucla.edu 
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The  central  schemes  we  use  in  this  paper  are  based  on  the  evolution  of  cell  averages  over  staggered  grids.  The 
staggered  versions  of  central  schemes  introduced  in  [19,  18,  11],  are  presented  in  §2.1  and  §2.2.  Central  schemes 
eliminate  the  need  for  a  detailed  knowledge  of  the  eigen-structure  of  the  Jacobian  matrices.  Instead  of  (ap¬ 
proximate)  Riemann  solvers  as  building  blocks  for  upwind  schemes,  simple  quadrature  formulae  are  used  for 
the  time  evolution  of  central  schemes.  This  approach  not  only  saves  the  costly  characteristic  decomposition  of 
the  Jacobians,  but  in  fact,  it  allows  us  to  completely  avoid  the  costly  evaluations  of  7  x  7  and  8x8  Jacobian 
matrices  in  one-  and  two-space  dimensions.  Moreover,  central  schemes  eliminate  the  need  for  dimensional  split¬ 
ting  which  is  particularly  relevant  for  the  multidimensional  MHD  system.  Indeed,  it  is  known  that  dimensional 
splitting  may  face  difficulties  with  the  propagation  of  other  than  genuinely  nonlinear  waves,  e.g.,  in  the  weakly 
hyperbolic  cases  reported  in  [13],  [8],  [11];  we  recall  that  the  MHD  equations  consist  of  such  waves  according 
to  a  main  observation  of  [2].  The  resulting  central  schemes  are  black-box,  Jacobian- free  MHD  solvers  whose 
sole  input  is  the  computed  MHD  fluxes.  The  fact  that  despite  their  simplicity  these  central  solvers  are  able  to 
resolve  accurately  the  complexity  of  one-  and  in  particular  two-dimensional  MHD  waves,  is  the  main  issue  of 
this  paper.  We  demonstrate  this  point  with  a  series  of  numerical  simulations. 

In  this  paper  we  focus  on  five  prototype  MHD  problems.  We  use  second-  and  third-order  non-oscillatory 
central  schemes  to  compute  the  approximate  solution  of  the  one-dimensional  Brio-Wu  shock  tube  models  [2]. 
Two  different  MHD  shock  tube  problems  are  studied  in  §3.1  and  §3.2.  The  second  order  Jiang-Tadmor  central 
scheme  [11]  is  implemented  for  the  approximate  solution  of  three  MHD  models  in  two  dimensions.  In  §4.1 
we  consider  the  Kelvin-Helmholtz  transverse  instability  problem  in  both  periodic  and  convective  formulations 
studied  earlier  by  e.g.,  [23,  12],  and  in  §4.2  we  study  the  Orszag-Tang  MHD  vortex  system,  [20],  [21],  Unlike 
the  Brio-Wu  and  Kelvin-Helmholtz  problems,  the  numerical  solution  of  the  Orszag-Tang  vortex  system  does 
not  necessarily  preserve  the  divergence  constraint,  V  •  B  =  0.  To  enforce  the  latter,  a  Leray  projection  corrector 
is  implemented  at  the  end  of  each  time-step,  replacing  the  computed  magnetic  field  with  its  divergence  free 
projection. 

Our  results  are  found  to  be  in  excellent  agreement  with  previous  simulations  of  the  same  problems  carried 
out  with  upwind-type  schemes,  [2,  23,  12],  and  complement  the  results  of  Wu  and  Chang  [25],  and  the  more 
recent  results  for  relativistic  MHD  flows  obtained  with  central-upwind  schemes  by  Del  Zanna  et.  al.  in  [6,  7]. 
These  results  demonstrate  the  ability  of  central  schemes  to  detect  and  resolve  the  discontinuous  solutions  that 
characterize  these  models,  while  retaining  efficiency  and  simplicity.  Indeed,  the  one-  and  two-dimensional 
reported  in  this  paper  demonstrate  the  remarkable  versatility  of  central  schemes  as  black-box,  Jacobian- free 
solvers  for  MHD  computations.  We  conclude  this  paper  with  a  two-page  appendix  which  provides  the  complete 
2D  central  code  for  the  convective  Kelvin-Helmholtz  problem. 

2  The  Numerical  Methods 

We  approximate  the  solution  of  (1.1) — (1.4)  using  predictor-corrector  central  schemes  which  are  implemented 
over  staggered  grids.  The  schemes  have  two  main  ingredients:  #1.  A  non-oscillatory  piecewise  polynomial 
reconstruction  of  pointvalues  from  their  cell  averages;  followed  by  #2.  Realizing  the  evolution  of  these  re¬ 
constructed  polynomials  in  terms  of  their  staggered  cell  averages.  Following  is  a  description  of  the  one-  and 
two-dimensional  schemes  together  with  the  magnetic  field  corrector  which  we  use  for  the  calculations  reported 
in  §3  (one-dimensional  results)  and  §4  (two-dimensional  results). 
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Figure  2.1:  Central  differencing  by  Godunov-type  scheme 

2.1  One  Dimensional  Methods 

We  begin  by  writing  the  system  of  equations  (1.1)-(1.4)  in  its  conservative  form 

ut  +  f(u)x  =  0,  (2.1) 

with 

u  =  (p,pvx,pvy,pvz,By,Bz,e)T ,  (2.2) 

f(u)  =  ( pvx ,  pv2x  +  p*  -  Bx,pvxvy  -  BxBy,  pvxvz  -  BXBZ,  Byvx  -  Bxvy, 

Bzvx  -  Bxvz,  (e  +  p*)vx  -  Bx{Bxvx  +  Byvy  +  Bzvz))T ,  (2.3) 

and  where  p*  =p  +  ^ B 2  stands  for  the  total  pressure  (static  plus  magnetic). 

Next,  we  follow  Godunov’s  seminal  idea  for  the  approximation  of  discontinuous  solutions  of  conservation  laws. 
Using  the  formulation  in  [18]  we  consider  the  sliding  averages,  u(x,t)  :=  Ix-Ax/2  u(£)  t)^£;  the  conservation 
law  (2.1)  then  reads 

w*(M)  +  ^:  f(u(x+  4r,t))  -f(u(x-  4r,t))  =o.  (2.4) 

Introducing  a  small  time  step  At,  and  integrating  over  the  slab  t  <  r  <  t  +  At  we  arrive  at 

I  r-t+M  a  a 

■u(x,  t  +  At)=  u( x,  t)  -  —  Jt  [f(u(x  +  — ,  t))  -  f(u(x  -  — ,  t))J  dr.  (2.5) 

So  far,  (2.5)  is  exact.  The  solution  to  (2.5)  is  now  realized  at  the  discrete  time  level  tn  =  nAt  by  a  piecewise 

polynomial  approximation,  w(x,tn)  ~  u(x,tn),  which  takes  the  form 

w(x,tn)  =  ^2pj(x)xj(x),  Xj{x)-=hr  (2.6) 

Here,  Pj(x)  are  algebraic  polynomials  supported  on  the  discrete  cells  Ij  =  IXj  =  [xj+i,Xj+i]  with  interfacing 
breakpoints  at  the  half-integers  gridpoints,  Xj±  i  =  (j  ±  |)Ax.  Sampling  (2.5)  at  x  =  Xj+i,  we  arrive  at  the 
new  staggered  cell  averages,  ,  centered  at  IJ+h=1^ 

1  1  r  1  [  rtn+1  rtTl+1 

=  -^ -J  w(x,tn)dx~—  J  f(w(xj+1,t))dt-  J  f(w(xj,t))dt  .  (2.7) 

The  evaluation  of  the  expressions  on  the  right  of  (2.7)  proceeds  in  two  steps,  which  will  occupy  the  rest  of  this 

section. 
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The  first  step  starts  with  the  known  cell  averages,  which  are  used  to  reconstruct  a  non-oscillatory 

piecewise  polynomial  approximation  w{x,tn)  =  Pj(x)xj{x).  We  use  piecewise  linear  functions  in  the  second 
order  case, 


Pj{x)  =  u^+w'jC- 


Ax 


l), 


and  piecewise  quadratic  functions  in  the  third  order  case, 


,,  „  //*^  7  \  1  ///*^  7  \  o 

Pj(x)  =  W3  +  WA~&r)  +  nwi  (  )  • 


Ax 


Here,  u>" ,  wr)  ,  and  w"  are  the  approximate  point  values  and  the  first  and  second  derivatives  of  w(x,tn)  at 
x  =  Xj,  which  are  reconstructed  from  the  given  cell  averages,  {w" }  7- .  Several  approximations  for  these  numerical 
derivatives  are  available  within  the  accuracy  constraints  of  the  schemes.  It  should  be  noted  that  the  procedure 
for  reconstruction  of  point  values  and  couple  of  numerical  derivatives  from  the  given  cell  averages  is  at  the 
heart  of  high-resolution,  non-oscillatory  central  schemes.  In  particular,  such  reconstructions  should  satisfy  the 
following  three  essential  properties: 

•  Vi  —  Conservation  of  cell  averages:  pj( x)\x—x.  =  w™ 

•  V2  Accuracy:  w(x,  tn)  =  u(x,  tn)+0((Ax)r )  for  r-  order  accurate  method,  wherever  u(-,t)  is  sufficiently 
smooth. 


•  V3  —  Non-oscillatory  behavior  of  JV  pj(x)xj(x)  which  is  characterized  in  different  ways  for  different 
reconstructions . 

Few  examples  are  in  order.  For  the  second-order  results  presented  in  §3,  we  use1 

w'j  =  MinMod(aA+iZt,-,  AoWj,  aA_Wj),  1  <  a  <  4.  (2-8) 

Here,  MinMod  stand  for  van-Leer’s  limiter,  [16],  where  MinMod(a,  6,  c)  =  sign(a)  min(|a|,  |6|,  |c|)  if  sign(a)  = 
sign{b)  =  sign(c),  and  it  vanishes  otherwise.  It  follows  that  this  reconstruction  procedure  is  non-oscillatory 
in  the  sense  of  satisfying  a  maximum  principle,  sup 3,  \  ^2jPj(x)Xj{x)\  <  supx  \^2jWj'Xj{x)\-  Moreover,  for  a 
restricted  set  of  a-values,  this  MinMod-based  reconstruction  is  Total- Variation  Diminishing  (TVD), 

II  Y^pAx)Xj{x)\\tv  <  II  Y1™?xAx)\\tv, 

3  3 

and  hence  the  corresponding  central  scheme  is  TVD,  [19].  For  the  third-order  case  we  use 


W3 

=  OjAow'j , 

(2.9) 

=  OjA+A-w™, 

(2.10) 

w " 

< 

=  w™ - L 

1  24’ 

(2.11) 

where  Oj  stands  for  the  third  order  limiters  specified  in  [18],  which  guarantee  the  non-oscillatory  behavior  of 
the  resulting  central  schemes,  this  time  with  a  Number  of  Extrema  Diminishing  (NED)  and  shape  preserving 
properties  (we  observe  that  starting  with  third-order  accuracy,  the  pointvalues  w™  differ  from  the  cell  averages, 
U!j  ).  Higher  order  non-oscillatory  extensions  are  offered  by  the  Essentially  non-oscillatory  (ENO)  reconstructions 
of  Harten  et.  al.,  e.g.,  [9],  and  their  implementation  with  central  schemes  can  be  found  in  e.g.,  [17]. 

Once  w(x,tn)  is  realized  as  a  piecewise  polynomial,  w{x,tn)  =  Y^Pj(x)Xj{x),  it  is  integrated  exactly  over  the 
interval  J3+i  to  compute  the  staggered  cell  averages  on  the  right  of  (2.7) 


1 

Ax 


w(x,  tn)dx 


1 

Ax 


pj+i(x)dx 


1  Where  A±  and  Aq  stand  for  the  usual  differences,  A ±Wj  =  ±(tCj±i  —  Wj ) ,  and  Aq 


^K  +  ^+1]  +  ^K-^+1]-(2.12) 

=  i(A+  +  A-) 
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We  turn  to  the  second  step  in  the  construction  of  the  central  scheme.  Here  we  follow  the  evolution  of  the 
reconstructed  point  values,  {w(xj,t  >  tn)}j,  along  the  midcells,  x  =  Xj,  which  are  governed  by 

wt  +  f(w)x  =  0,  t  >  tn\  w(x,tn)  =  pj(x),x  €  Ij  (2-13) 

If  {a,k(u)}k  are  the  eigenvalues  of  the  Jacobian  A(u)  :=  then  by  hyperbolicity,  information  regarding  the 
interfacing  discontinuities  at  {xj±  1 ,  tn)  propagates  no  faster  than  maxj  |ofc(u)|.  Hence,  the  midcell  values, 
{w(xj,r  >  tn)}j ,  remain  free  of  discontinuities  as  long  as  the  CFL  condition,  A t  <  \  -  maxj.  |afc(«)|,  is  met. 
Therefore,  the  flux  integrals  in  the  right  of  (2.7)  involve  only  smooth  integrands  and  can  be  evaluated  with 
appropriate  quadrature  rules  to  any  desired  degree  of  accuracy.  In  particular,  the  second  order  Nessyahu- 
Tadmor  scheme  [19]  makes  use  of  the  midpoint  rule 


f  [w  (xj ,  t)  )  dr 


A  i/0 


n+i 


), 


while  the  third  order  Liu-Tadmor  scheme  [18]  makes  use  of  Simpson’s  rule 


f(w{Xj,T))dT 


At 


/(<)  +  4/(<H’)  +  /(<+1) 


These  quadrature  formulae  require  the  computation  of  the  intermediate  point  values  w"+/3,  /3  =  0,  1.  A 

natural  approach  for  computing  these  point  values  employs  Taylor’s  expansion  and  the  differential  equation 
(2.1),  wt  =  —f(w)x.  The  resulting  predictor  step  in  the  second  order  case  reads 


wTl  =  U  =  /«)■  (2-14) 

Here  A  =  At/ Ax  is  the  fixed  mesh  ratio  and  /'  is  an  approximate  numerical  derivative  —  an  approximation 
to  the  spatial  derivative  ft/ Ax  ss  f(wf)x.  Different  recipes  for  numerical  derivatives  are  available  within  the 
accuracy  constraints  of  our  calculations.  We  shall  mention  two  of  them  which  retain  the  overall  second-order 
accuracy,  consult  [19].  In  the  first  approach  we  utilize  the  chain  rule,  f{w)x  =  A(w)wx.  With  A{w ™)  standing 
for  the  computed  Jacobian  and  w/  denoting  the  reconstructed  numerical  derivative,  vj/ /Ax  «  w(xj ,tn)x,  we 
have 


fj  =  A(w’j)w'j,  w'j  =  MinMod(aA+u>j,  AoWj,  aA_Wj).  (2-15) 

Alternatively,  we  can  implement  a  spatial  MinMod  limiter  directly  for  numerical  differentiation  of  the  grid 
function  {fj  =  f(wf)}j, 

fj  =  MinMod(aA+/j,  A0/j,  aA_/j),  fj  =  /«).  (2.16) 

Remark  that  the  resulting  predictor-step  (2. 14), (2. 16)  completely  avoids  the  evaluation  of  the  Jacobian  A(w "). 
Similar  recipes  are  available  for  the  third  order  predictor, 

WT0  =  WJ  -  W(w7  Yf,i)y’  (2-17) 

As  before,  the  evaluation  of  the  expression  on  the  right  of  (2.17)  can  utilize  the  chain  rule, 


f(w  -  =  f(w  ~  ~YA^Wx\x  =  A(w)wx  ~  ~Y (^2(w)^x)x  +  0(X)2. 


A/3 


A/3 


Using  the  third-order  accurate  numerical  derivative  in  (2.9),  we  end  up  with  two  possible  recipes  for  computing 

—  f') 

J  2  ■’]> 


the  third-order  accurate  predicted  mid-values,  { f(w] /  —  in  (2.17), 


{/K  -  ¥ffj)Y  =  0jAo {/K-^^KK)},  /?  =  w^OjAoWj, 
{/«  -  -  A(wj)w'j  -  ?Y9jAo{A(wj)2Wj},  /3=  |,1,  w/  =  e3A0w3. 
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(2.18) 

(2.19) 


A  third  possibility,  a  so-called  Jacobian- free-form  (JFF)  alternative  is  provided  by 

{/K-^/i)}'  =  ^Ao(/(<-^A0(/K)))),  /?=±1.  (2.20) 

Again,  the  computation  of  the  intermediate  values  w™+l3  in  (2.20)  avoids  the  explicit  evaluation  of  the 
Jacobians,  A(u>”).  The  high-resolution  of  these  Jacobian-free  versions  (2.16)  and  (2.20)  is  evident  in  the 
numerical  results  reported  in  §3.  The  intermediate  values  can  also  be  approximated  by  Runge-Kutta 

solvers  of  the  ODE  wr  =  fx\x=Xj,  w(xj,0)  =  w™,  r  >  tn,  where  fx  stands  for  the  numerical  derivative  of  /.  This 
approach  offers  yet  another  avenue  for  the  time  evolution  of  central  schemes:  in  particular,  we  mention  in  this 
context  the  efficient  evolution  procedure  based  on  higher  order  Runge-Kutta  Natural  Continuous  Extensions, 
e.g.,  [3,  17],  which  significantly  reduces  the  number  of  computations  per  grid  point. 

Equipped  with  the  predicted  point  values,  we  are  set  to  evaluate  the  integrals  on  the  right  hand  side  of  (2.7). 
They  are  approximated  by  the  midpoint  rule  in  the  second  order  case, 

Ai/(«£+*)  =:  Af/;+5,  (2.21) 

and  by  Simpson’s  rule  in  the  third  order  case, 

rtn+1  a  t  i  i 

j  f{w{xj,T))dT  «  — [/«)  +  4/>;+5)  +  fiwp1)]  =:  A tf;+\  (2.22) 

n+i 

Denoting  these  approximate  values  of  the  flux  integrals  by  /■  2 ,  the  corrector  step  for  both  the  second-  and 

third-order  central  schemes  amount  to  the  same  statement, 


f(w(Xj,T))dT  : 


uj+ 1 


3+ 1 


-A 


„n+i  j .n+  = 

fj+ 1  -  fj 


(2.23) 


2.2  Two  Dimensional  Method 

We  begin  the  description  of  the  two-dimensional  scheme  by  writing  (1.1)-(1.4)  in  conservation  form 

ut  +  f(u)x  +  g(u)z  =  0,  (2.24) 

where 

u  =  (p,pvx,pvy,pvz,Bx,By,B,,e)T,  (2.25) 

f(u)  =  (pvx,pvl  +  p*  -  B2x,  pvxvy  -  BxBy ,  pvxvz  -  BXBZ,  0,  Byvx  -  Bxvy, 

Bzvx  -  Bxvz,  (e  +  p*)vx  -  Bx(Bxvx  +  Byvy  +  Bzvz))T ,  (2.26) 

!J ( c)  (pc2 .  p‘OzOx  BzBx.  pvzvy  BzBy,  pvz  T  p  Bz1Bxvz  B z vx , 

Byvz  -  Bzvv,  0,  (e  +  p*)vz  -  Bz(Bxvx  +  Byvy  +  Bzvz))T .  (2.27) 

We  introduce  space  scales  Ax  and  A z,  and  consider  the  sliding  averages  of  (2.24)  over  the  two  dimensional  cell 
[x  -  ^,x  +  A*]  x  [z  -  As,  2  +  As] 

1  rz~\  2~  /\  71 

ut{x,z,t)  +  [/(u(a;+ —  ,??,t))  - /(w(a;  -  —  ,r?,t))]d?7 

"'2  T 

+  Ax  A-  /  »t))]^  =  °  (2-28) 

We  proceed  along  the  same  lines  of  the  one  dimensional  case;  integration  over  the  slab  t  <  r  <  t  +  At  is 
performed,  and  the  solution  of  the  resulting  equation  is  approximated  at  time  level  tn  —  nAt  by  a  piecewise 
bilinear  polynomial  w(x,z,tn )  =  Y^Pjk(x,  z)xjk(x,  z),  where  Xjk{x,z )  =  1/  x/*fc  is  the  characteristic  function 
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of  the  control  volume  IXj  x  Izk  •  Choosing  x  =  Xj+ 1/2,  and  2  =  zk+i/2  as  the  interfacing  boundaries,  the  two 
dimensional  version  of  equation  (2.7)  reads 


Wj+i,k+i  wj+i,k+  i 


1 

AxAz 

1 

AxAz 


xj+ 1 


[^(w(a;,  2fc+i,t))  -  g(w(x,  2/s,i))]cfadt. 


(2.29) 


The  reconstruction  of  the  point  values  w(x,z,tn)  is  carried  out  using  piecewise  bilinear  functions  pjk(x,z)  = 
Wj'k  +  u>jk  ( X ^x3 )  +  Wjk  ( ) ,  with  the  approximate  numerical  derivatives  in  the  x  and  2  directions  given  by 


w'jk  =  MinMod(aA+xw;"fc ,  A0xwJk ,  aA_xWjk), 
w)k  =  MmMod(aA+zWjk,A0zWjk,aA_zw™k). 

These  piecewise  polynomials  admit  the  following  cell  averages  on  the  right  of  (2.29) 


rxi+  i  Czk+z 


wj+i,k+i 


Pjk{x,  z)dzdx 


rxi+ 1  rzk+  * 


Pj+i,k{x,  z)dzdx 


'X.,1  J  Zk 


rxi+ i  /***+! 


Pj,fc+i(a;,  z)dzdx  + 


rxj+i  rZk+l 


Pj+i,k+i(x,  z)dzdx 


'Xj+h  JZk+h 


—  ~^(W7k  +  W7+l,k  +  W7,k+ 1  +  ^j+l.fc+l) 


16  L 


(w'fc  -  w'+1,fc)  +  (^,fc+l  -  w'  +  l,fc+l)  +  (w}fe  -  W},fc+l)  +  (w)+l,k  -  w)+ l,fc+l) 


(2.30) 

(2.31) 


(2.32) 


The  reconstruction  procedure  described  above  enjoys  the  conservation,  accuracy  and  non-oscillatory  properties, 
V1-V3,  analogous  to  those  listed  for  one  dimensional  schemes  in  §2.1.  In  particular,  the  non-oscillatory  property 
in  this  case  is  characterized  by  the  scalar  maximum  principle,  consult  [11]. 

To  approximate  the  flux  integrals  on  the  right  of  (2.29),  we  first  predict  the  midpoint  values 


n+i 

wit 


(2.33) 


where  A  =  At/ Ax  and  p  =  At/ Az  are  the  (fixed)  mesh  ratios  in  the  x  and  2  directions,  and  f/k  and  y)fe  stand 
for  the  numerical  slopes  of  f(w)  and  g(w).  In  this  case,  we  choose  the  Jacobian  free  form 


f'jk  =  MinMod(aA+xfjk,A()xfjk,aA_xfjk),  (2.34) 

g)k  =  MinMod(aA+z0j-fc,  A0zgjk,aA_zgjk).  (2.35) 


In  the  computations  below  we  set  the  free  parameter  a  =  1.4.  The  flux  integrals  are  then  approximated  by  the 
midpoint  rule  in  time,  and  by  second  order  rectangular  quadrature  in  space 


Zfc+l 


f(w(xj+i,z,  t))dzdt < 


Xj  + 1 


g(w(x,  Zk+i,  t))dxdt  • 


AzAt 


f(wjti,k)  +  /(w”+i%+i) 


Ax  At 


9  (' 


n+i 

yi+i,fc+i 


)  +  9( 


n+h  \ 

u,-,fc+i) 


(2.36) 

(2.37) 
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The  new  staggered  cell  averages  k+i  }j,fc  read 


1 

16 

1 


j+l,fc  '  wj,k-\- 1 

A  r  .  n+  J 

—  ft  in 

2  P 


'  ^i+ijfe+i  J 
n+J 


—  (w'fc  -  ic'+1>fc)  -  -  /(u>j+1%)  -  f(wjk  2) 


'j+i,fc+i)  2  L 


/(™"+l%+l)  -  /(«£fc+l) 


n-H 


iK-fe  -  <*,+1)  -  t  fsKfe)  -  5(^”fc+5) 


16 v  J'fc 

1 


2  L 


■  -*■  /  \  \  \  A*  /  «-i-f  \  /  «-i-f  \ 

+  -r^K+i.fc  -  ^j+i,fc+i)  -  9  9{Wj+1,k+i)  -  g{wj+ iifc) 


n+J 


n-H 


16 


2  L 


(2.38) 


In  addition  to  reconstruction  and  evolution,  our  second-order  two-dimensional  scheme  requires  the  correction 
of  the  magnetic  field,  B,  in  order  to  guarantee  the  solenoidal  constraint  V  •  B  =0.  We  satisfy  this  constraint 
by  updating  the  cell  averages  of  the  magnetic  field  at  the  end  of  each  time  step,  replacing  the  computed  B, 
with  its  divergence- free  projection,  Be.  To  this  end  the  so  called  Leray  projection  is  carried  out  by  solving  the 
Poisson  equation 


A^=-V-B  (2.39) 

with  the  appropriate  boundary  conditions.  Since  the  coplanar  structure  of  the  problem  guarantees  dBv/dy  =  0, 
only  the  components  Bx  and  Bz  need  to  be  updated.  We  use  a  fast  Poisson  solver  for  the  standard  five  point 
discretization  of  the  potential  (f>,  and  central  differences  for  the  divergence  of  the  magnetic  field,  V  •  B.  The 
corrected-divergence  free-magnetic  field,  Bc,  is  then  realized  as 

Bc  =  B  +  V<jf>.  (2.40) 

Applying  the  divergence  operator  V  =  (d/dx,  d/dz)  to  both  sides  of  (2.40),  one  can  easily  verify  that  V-Bc  =  0. 
We  use  here  only  one  out  of  several  methods  to  enforce  the  so-called  ‘constrained  transport’  and  we  refer  the 
reader  to  [1],  [22]  and  the  references  therein  for  a  general  discussion  and  [10]  for  handling  the  solenoidality 
constraint  in  the  context  of  MHD  schemes  over  staggered  grids. 

3  One  Dimensional  Numerical  Results 

In  this  section  we  present  numerical  simulations  of  the  one-dimensional  MHD  equations  (2.1)-(2.3).  The  results 
were  obtained  using  different  versions  of  the  second-order  Nessyahu-Tadmor  central  scheme  (2.14)-(2.16),  (2.21), 
(2.23),  based  on  the  minmod  limiter  (2.8),  and  third-order  Liu-Tadmor  central  scheme,  (2.17)-(2.20),  (2.22)- 
(2.23),  based  on  the  non-oscillatory  limiter  (2.9)-(2.11).  The  schemes  are  implemented  for  the  approximate 
solution  of  two  coplanar  shock  tube  MHD  models  described  by  Brio  and  Wu  in  [2].  We  use  a  uniform  grid  in 
the  space  discretization,  and  in  both  cases  we  choose  the  time  step  dynamically  with  CFL  restriction 


At  =  0.4 


Ax 

maxfc  |ofe(u)| 


where  {a*, (it)} a-  are  the  eigenvalues  of  the  Jacobian  matrix  of  /(it). 


(3.1) 


3.1  Brio— Wu  Shock  Tube  Problem 


The  first  one  dimensional  Riemann  problem  we  consider  consists  of  a  shock  tube  with  two  initial  equilibrium 
states,  ui  and  ur, 


{pi  Vx  lVyiVz,  By  ,  Bz ,  p) 


(1.0,  0,0,0, 1.0, 0,1.0)T  for  a;  <  0, 
(0.125,  0,  0, 0,  -1.0,  0, 0.1)T  for  x  >  0, 


(3.2) 


and  complemented  with  the  constant  values  of  Bx  =  0.75  and  7  =  2.  The  problem  is  solved  for  x  £  [— 1, 1] 
with  800  grid  points,  and  numerical  results  are  presented  at  t  =  0.2.  Figures  3. 2-3. 5  show  the  density,  the 


x —  and  y— velocity  components,  the  y— magnetic  field,  and  pressure  profiles  computed  with  different  choices  of 
numerical  derivatives  /j  and  u'j. 

We  note  that  the  hydro  dynamical  data  of  this  problem  is  the  same  as  that  in  Sod’s  shock  tube  problem  of 
gas  dynamics.  The  variety  of  MHD  waves,  however,  poses  a  considerable  challenge  for  high-resolution  methods 
such  as  the  JFF  central  schemes  described  in  this  paper.  The  solution  of  this  problem  consists  of  a  left-moving 
fast  rarefaction  wave  (FR),  a  slow  compound  wave  (SM)  which  results  from  an  intermediate  shock  that  changes 
By  from  0.58  to  —0.31  and  a  slow  rarefaction  that  changes  By  from  —0.31  —0.53,  a  contact  discontinuity  (C), 
a  right-moving  slow  shock  (SS),  and  a  right-moving  fast  rarefaction  wave  (FR).  Note  that  the  solution  to  this 
problem  is  not  unique  if  Bz  and  vz  are  not  identically  zero. 

Figures  3. 2-3. 5  show  the  solutions  calculated  with  different  versions  of  the  second  and  third  order  schemes. 
These  results  are  comparable  with  the  second  order  upwind  computations  of  Brio  and  Wu  in  [2] ,  and  with  the 
fifth  order  WENO  computations  presented  by  Jiang  and  Wu  in  [12].  Our  numerical  results  demonstrate  that 
central  schemes  —  while  avoiding  any  characteristic  information  other  than  an  estimate  of  the  maximal  speed 
max  a-  |afc(u)|,  they  still  capture  the  main  features  of  the  discontinuous  MHD  solutions. 

In  figures  3.4  and  3.5  we  observe  the  third  order  oscillations  near  the  trailing  edge  of  the  fast  rarefaction  wave 
that  are  less  evident  with  the  second  order  results.  This  is  due  to  the  higher  order  polynomial  reconstruction; 
indeed,  the  same  phenomenon  is  reported  in  [12]  when  comparing  fifth  order  WENO  with  second  order  results.  A 
final  remark  regarding  the  Jacobian  Free  Forms  (JFF),  (2.16),  (2.20):  not  only  that  they  offer  a  more  economical 
approach  by  avoiding  costly  matrix  multiplication,  but  they  also  provide  a  better  control  of  these  oscillations 
as  well  as  a  better  resolution  of  the  contact  discontinuity  (C);  in  particular,  they  allow  the  expected  jump  in 
the  density  profile  without  disrupting  the  constant  state  of  the  remaining  conserved  quantities. 


Figure  3.2:  Results  of  Brio-Wu  shock  tube  problem  at  t  =  0.2  computed  with  800  grid  points  using  the  second 
order  central  scheme.  Here  /(  is  computed  componentwise  with  the  minmod  limiter  using  the  JFF  version 
(2.16)  with  a  =  1.4. 
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Figure  3.3:  Results  of  Brio-Wu  shock  tube  problem  at  t  =  0.2  computed  with  800  grid  points  using  the  second 
order  central  scheme.  Here,  the  numerical  flux,  /j  =  A(w^)w'y  is  evaluated  using  the  minmod  limiter  (2.15) 
with  a  =  1.4. 


Figure  3.4:  Results  of  Brio-Wu  shock  tube  problem  at  t  =  0.2  computed  with  800  grid  points  using  the  third 
order  central  scheme.  Here  /j  are  computed  component-wise  using  the  JFF  version  (2.20)  with  a  =  1.4. 
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Figure  3.5:  Results  of  Brio-Wu  shock  tube  problem  at  t  =  0.2  computed  with  800  grid  points  using  the  third 
order  central  scheme.  Here,  the  intermediate  numerical  fluxes,  {f(wj  —  A/3/2 fj)}',  f3  =  0,  |,1,  are  evaluated 
using  the  third-order  limiter  (2.18)  with  a  =  1.4. 


3.2  Brio— Wu  High  Mach  Shock  Tube  Problem 

The  following  shock  tube  model  proposed  by  Brio  and  Wu  in  [2],  is  used  to  check  the  robustness  of  the  numerical 
schemes  for  high  Mach  number  problems.  The  initial  equilibrium  states,  ui  and  ur,  are  given  by 


I  (1.0, 0,0,0, 1.0,  0,1000)T  for  £  <  0, 
\  (0.125, 0,0,0, -1.0, 0,0. 1)T  for  £  >  0, 


(3.3) 


complemented  with  the  values  of  Bx  =  0,  and  7  =  2.  The  Mach  number  of  the  right-moving  shock  wave  is 
15.5.  If  the  plasma  pressure  is  replaced  by  the  sum  of  the  static  and  magnetic  pressures-denoted  by  p*  above, 
the  problem  becomes  a  standard  hydrodynamical  Riemann  problem.  The  solution  is  presented  at  t  =  .012, 
x  £  [—1, 1],  with  200  grid  points  and  with  CFL  number  0.4,  consult  (3.1). 
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Figure  3.6:  Results  of  Brio-Wu  high  Mach  problem  at  t  =  0.012  computed  with  200  grid  points  using  second 
order  central  scheme.  Here  /j  are  computed  component-wise  using  the  JFF  (2.16)  with  a  =  1.4. 


Figure  3.7:  Results  of  Brio-Wu  high  Mach  problem  at  t  =  0.012  computed  with  200  grid  points  using  third 
order  central  scheme.  Here,  the  intermediate  numerical  fluxes,  {/('«>,  —  A/3/2/j)}',  (3  =  0,|,1,  are  evaluated 
using  the  third-order  limiter  (2.18)  with  a  —  1.4. 


12 


The  solution  of  this  second  Riemann  problem  consists  of  a  left-moving  fast  rarefaction  wave  (FR),  followed 
by  a  tangential  discontinuity  (TD),  and  a  right  moving  fast  shock  (FS)  with  Mach  number  15.5.  Across  the 
tangential  discontinuity,  the  density,  the  magnetic  field  and  the  pressure  can  change,  but  both,  the  fluid  velocity 
and  the  total  pressure  are  continuous. 

As  in  the  previous  problem,  our  results  are  comparable  to  those  obtained  by  Brio  and  Wu  in  [2]  with  their 
second  order  upwind  method  and  the  ones  presented  by  Jiang  and  Wu  in  [12],  computed  with  a  fifth  order 
WENO  scheme.  These  results  indicate  that  the  methods  described  above  are  robust  even  in  their  ‘greedy’, 
Jacobian-free  version. 


4  Two  Dimensional  Numerical  Results 

In  this  section  we  present  the  numerical  solutions  of  two  prototype  models  of  two-dimensional  MHD  equations. 
For  the  first  problem  —  the  Kelvin-Helmholtz  instability  with  transverse  magnetic  field  configuration,  we 
consider  two  different  sets  of  boundary  conditions  in  the  direction:  periodic  in  the  first  case  and  a  free 
outflow  boundary  in  the  second  convective  setup.  The  second  Orszag-Tang  problem  introduced  by  Orszag 
and  Tang  in  [20]  as  a  simple  model  to  study  MHD  turbulence,  and  has  become  a  standard  model  to  validate 
numerical  algorithms.  In  both  cases,  the  time  scale,  At  is  determined  dynamically  with  CFL=0.4, 


At  = 


0.4 

(maxfc  |ak(«)|/Aar)  +  (maxfc  \bk(u)\/Az) 


(4.1) 


where  {ofc(u)}fc  and  {6fc(u)}fc  represent  the  eigenvalues  of  the  Jacobian  matrices  of  f(u)  and  g(u)  respectively. 


4.1  Transverse  Kelvin-Helmholtz  Instability 

The  Kelvin-Helmholtz  instability  arises  when  two  superposed  fluids  flow  one  over  the  other  with  a  relative  veloc¬ 
ity.  It  models,  for  example,  the  important  mechanism  for  the  momentum  transfer  at  the  Earth’s  magnetopause 
boundary,  which  separates  the  solar  wind  from  the  Earth’s  magnetosphere,  [12].  We  apply  the  second  order 
central  scheme  of  Jiang  and  Tadmor,  (2.33)-(2.38),  for  the  two  dimensional  periodic  and  convective  models 
with  transverse  magnetic  field  configuration. 

In  both  cases,  the  governing  equations  are  (2.24)-(2.27)  are  subject  to  initial  conditions 

(p,  vx,  vy,  vz,Bx,  By,  Bz,p)T  =  (1.0,  vxo  +  vx0, 0, 0, 0, 1.0, 0, 0.5)T,  (4.2) 


where 


VxO 


VxO 


and 


f  -v0  sin 

l  0, 


if  —  #  <  ®  <  # 

otherwise, 


(4.3) 

(4.4) 


with  Vo  =  2,  Vq  =  —0.008,  A  =  57r  and  a  =  1.  Also,  the  grids  are  stretched  in  the  z— direction  with  a  Roberts 
transformation 


H  sinh  (tz/2H) 
sinh  (r/2) 


t=  6, 


(4.5) 


which  renders  a  denser  grid  near  z  =  0  where  the  effect  of  the  small  initial  perturbation  vxo  is  more  noticeable, 
and  a  coarser  grid  near  z  =  ±U,  where  little  action  takes  place. 

In  the  periodic  case,  the  computational  domain  is  [— ,  4]  x  [0:  H\,  with  L  =  57r  and  H  =  10.  The  free  outflow 
condition  is  applied  at  the  top  boundary,  z  =  H ,  and  the  bottom  boundary  values  are  recovered  by  symmetry, 
since  p,  p ,  and  By  are  symmetric,  and  vx  and  vz  are  antisymmetric  under  the  transformation  (x,  z)  — >  (— x,  —z). 
In  figure  (4.8),  we  present  solutions  at  t  =  144,  with  96  x  96  and  192  x  192  grid  points. 

The  resolution  and  accuracy  of  our  results  are  comparable  to  those  achieved  by  the  second  order  gas  kinetic 
scheme  of  Tang  and  Xu,  [23]  using  a  200  x  200  uniform  grid,  and  to  those  achieved  by  Jiang  and  Wu’s  fifth 
order  WENO  scheme,  [12],  obtained  with  a  coarser  grid.  It  should  be  noted,  however,  that  in  the  case  of  central 
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schemes,  the  additional  computational  cost  generated  by  thinner  grids  is  compensated  by  the  simplicity  of  the 
algorithm:  no  characteristic  decompositions  are  computed,  no  Jacobians  are  required  and  dimensional  splitting 
is  avoided. 

In  the  convective  setup,  the  initial  conditions  and  perturbation  are  the  same  as  in  the  periodic  setup,  (4.2)- 
(4.4).  In  this  case,  the  free  outflow  condition  is  applied  to  all  four  boundaries  of  the  computational  domain 
[—  -j,  j\  x  [—  H1  H]\  Here  H  =  20  and  L  =  557r,  with  L  »  A-so  chosen  to  allow  the  excitation  to  convect  freely 
without  disturbing  the  x— boundaries. 

Figures  4.9  and  4.10  display  the  solution  computed  with  1056  x  192  grid  points  at  t  =  120  and  t  =  145 
respectively.  As  in  the  periodic  case,  our  results  are  comparable  to  those  previously  obtained  with  upwind 
schemes,  [12,  23],  and  demonstrate  the  ability  of  central  schemes  to  detect  and  resolve  steep  gradients  without 
any  detailed  knowledge  of  the  eigen  structure  of  the  Jacobian  matrix. 


Figure  4.8:  Results  of  transverse  Kelvin-Helmholtz  instability  with  periodic  x— boundary  conditions.  Left 
column  uses  96  x  96  points  and  right  column  uses  192  x  192  grids  respectively.  There  are  20  contours  for  density 
and  pressure.  Red-high  value,  blue-low  value.  Density  ranges  from  0.79  to  1.2,  pressure  range  from  0.32  to 
0.71  and  maximum  value  for  the  velocity  is  1.25. 
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Figure  4.9:  Solution  of  convective  Kelvin-Helmholtz  instability  at  t=120  with  1056  x  192  grid  points.  The 
density  ranges  from  0.63  to  1.3,  and  the  pressure  ranges  from  0.20  to  0.85,  the  maximum  value  for  the  velocity 
is  1.54. 
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Figure  4.10:  Solution  of  convective  Kelvin-Helmholtz  instability  at  t=145,  with  1056  x  192  grid  points.  The 
Density  ranges  from  0.43  to  1.3,  pressure  ranges  from  0.10  to  0.86  and  maximum  value  for  the  velocity  is  1.94. 


Figure  4.11:  Time  evolution  of  the  total  transverse  kinetic  energy,  log(^  f  pv^dxdz),  integrated  over 
[—  L/2,L/2]  x  for  both  periodic  and  convective  Kelvin-Helmholtz  instability.  The  results  for  the 

periodic  case  with  96  x  96  and  192  x  192  grid  points  are  represented  by  a  dashed  and  a  doted  curve  respectively. 
The  convective  configuration  is  represented  by  a  solid  line. 
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4.2  Orszag-Tang  MHD  Turbulence  Problem 

This  model  considers  the  evolution  of  a  compressible  Orszag-Tang  vortex  system.  The  evolution  of  the  vortex 
system  involves  the  interaction  between  several  shock  waves  traveling  at  various  speed  regimes  [12,  24],  which 
makes  the  problem  especially  attractive  for  numerical  experiments. 

The  initial  data  is  given  by 

p(x,  z,  0)  =  y2,  vx(x,  z,  0)  =  —  sin  z,  vz(x,  z,  0)  =  sin  x, 
p(x,z,  0)  =  7,  Bx{x,  z,  0)  =  —  sin  z,  Bz(x,  z,  0)  =  sin  2x 

where  7  =  5/3.  With  this  initial  data,  the  root  mean  square  values  of  the  velocity  and  magnetic  fields  are  both 
1;  the  initial  average  Mach  number  is  1,  and  the  average  plasma  beta  is  10/3. 

We  solve  the  problem  in  [0,  27t]  x  [0,  27t],  with  periodic  boundary  conditions  in  both  x-  and  z-directions  using 
a  uniform  grid  with  384  x  384  points. 

In  this  problem,  the  numerical  scheme  fails  to  satisfy  the  divergence  free  constraint  of  the  magnetic  field, 
V  •  B  =  0.  In  order  to  guarantee  this  condition  and  avoid  numerical  instability,  we  project  the  updated  magnetic 
field  B  into  its  divergence  free  component  at  the  end  of  every  time  iteration  by  applying  the  correction  (2.39)- 
(2.40). 


Figure  4.12:  Orszag-Tang  MHD  turbulence  problem  with  a  384  x  384  uniform  grid  at  t  =  0.5.  There  are  12 
contours  for  density  and  pressure.  Red-high  value,  blue-low  value.  Density  range  from  2.1  to  5.8,  pressure 
range  from  1.0  to  5.7.  The  maximum  values  of  |r>|  and  \B\  are  1.6  and  1.6  respectively. 

Figures  (4.12),  (4.13),  and  (4.14)  present  the  solution  of  the  Orszag-Tang  vortex  system  at  t  =  0.5,  t  =  2, 
and  t  =  3  respectively.  Again,  these  results  are  comparable  with  those  obtained  by  upwind  schemes.  The 
lower  order  of  our  methods  require  a  thinner  grid  than  Jiang  and  Wu’s  fifth  order  WENO  scheme  for  resolving 
accurately  the  various  shocks  that  the  vortex  system  develops.  However,  this  relative  loss  of  efficiency  is  still 
compensated  by  the  simplicity  in  the  implementation  of  the  schemes. 
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Figure  4.13:  Orszag-Tang  MHD  turbulence  problem  with  a  384  x  384  uniform  grid  at  t  =  2.  There  are  12 
contours  for  density  and  pressure.  Red-high  value,  blue-low  value.  Density  ranges  from  0.62  to  6.3,  pressure 
ranges  from  0.14  to  7.0.  The  maximum  values  of  |u|  and  \B\  are  1.6  and  2.8  respectively. 


Figure  4.14:  Orszag-Tang  MHD  turbulence  problem  with  a  384  x  384  uniform  grid  at  t  =  3.  There  are  12 
contours  for  density  and  pressure.  Red-high  value,  blue-low  value.  Density  ranges  from  1.2  to  6.1,  pressure 
ranges  from  0.34  to  6.3.  The  maximum  values  of  |r>|  and  \B\  are  1.7  and  3.0  respectively. 
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Figure  4.15:  Pressure  distribution  along  the  line  z  =  0.6257T  at  t  =  3. 


5  Conclusions 

The  numerical  results  presented  in  this  paper  demonstrate  the  ability  of  central  schemes  to  compute  the  dis¬ 
continuous  solutions  of  ideal  MHD  equations  accurately.  Our  numerical  tests  are  in  excellent  agreement  with 
the  results  of  Jiang  and  Wu  [12],  and  they  complement  previous  results  of  Wu-Chang  [25]  and  the  more  recent 
results  of  Del  Zanna  et.  al.  [6,  7]  where  the  efficiency  of  central  schemes  is  demonstrated  for  other  MHD  models. 

How  one  can  quantify  this  efficiency?  the  answer  depends  on  too  many  local  features  —  computer  code, 
hardware  and  database  configurations,  ...  which  prevent  a  precise  quantitative  answer.  We  therefore  report 
here  on  our  ‘subjective’  results  of  CPU  running  time  for  solving  the  above  Orszag-Tang  problem,  comparing  the 
second-order  staggered  central  scheme  [11]  vs.  WENO  scheme  [12]  (fifth  order  in  space  and  forth  in  time).  The 
Jacobian-free  form  (JFF)  of  the  second-order  central  schemes  offered  a  speed-up  factor  of  25  in  this  case.  We 
should  point  out,  however,  that  the  second-order  results  required  a  refinement  of  the  spatial  grid  by  a  factor  of 
~  2,  in  order  to  achieve  a  resolution  similar  to  the  fifth  order  WENO.  With  the  more  restrictive  CFL  condition 

the  staggered  central  CFL  condition  is  ~  0.5  rather  than  1,  this  increases  the  amount  of  2D  ‘work’  by  a  factor 
8.  Thus,  the  acceleration  factor  for  a  given  resolution  in  this  case  is  a  factor  of  ~  3.  The  additional  gain  offered 
by  the  black-box  central  solvers  lies  in  their  simplicity:  neither  characteristic  decomposition,  nor  dimensional 
splitting  is  required.  The  relative  ease  of  implementation  is  highlighted  by  the  2D  code  in  the  appendix  below, 
where  the  intricate  eigensystem  specified  in  [12,  pp.  570-572]  is  completely  avoided. 

Finally,  we  observe  that  the  lower  resolution  of  our  staggered  central  schemes  necessitates  further  refinement 
of  the  spatial  grids  in  order  to  achieve  the  same  resolution  as  higher  order  upwind  methods,  resulting  in  a  loss  of 
efficiency.  We  will  address  this  issue  in  a  second  part  of  our  work  on  central  schemes  for  MHD  equations,  where 
we  use  a  semi-discrete  version  of  the  central  schemes  introduced  in  [14,  15]  and  its  higher-order  extensions, 
consult  [3,  17].  This  approach  retains  both  the  simplicity  of  implementation  and  the  efficiency  of  Riemann- 
solver-free  algorithms. 


A  Appendix.  A  Central  Code  for  Two-Dimensional  Ideal  MHD 
Equations 

The  following  C++  code  was  used  to  compute  the  solution  of  the  convective  Kelvin-Helmholtz  Instability,  figures 
(4.9)-(4.11).  We  refer  the  reader  to  http://www.math.ucla.edu/~jbalbas/MHD/code  for  this  two  dimensional 
MHD  code  and  additional  software  required  for  its  implementation. 

///////////////////////////////////////////////////////////////  //  type  reconstruction 

//  2nd  order  Central  Scheme  for  2-d  MHD  Equations:  // 

//  //  Input:  grid  size  J  x  K,  limiter  parameter  alpha,  CFL  number. 

//  Transverse  Kelvin-Helmholtz  instability — convective  set-up.  // 

//  //  Output:  density,  velocity  field,  magnetic  field,  pressure. 

//  Jacobian  Free  form  of  Staggered  LxF  with  MUSCL  /////////////////////////////////////////////////////////////// 
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CAMdoubleVector  e_int (42000) ,  time (42000); 


int  mainO 

////////////////////////////////////////////////////////// 

//  Grid 

////////////////////////////////////////////////////////// 
double  x_init,  x_final,  z_init,  z_final,  t_init,  t_final, 
double  dx,  dz_eq,  dt,  lambda; 
double  alpha,  pi; 
double  gamma=2.0; 
double  cfl; 

double  rho_doub,  p_doub,  A,  B,  max_eigx,  max_eigz; 
double  e_trans; 
long  J,  K; 

long  j ,  k,  jl,  j2 ,  1; 
long  mid,  s=0; 
bool  odd=true; 

if stream  InFile; 

InFile . open ("input " ,  ios::in); 
if  (! InFile) 

cerr«"Error  in  opening  input  file"; 
exit (1)  ; 

> 

InFile»J ; 

InFile»K; 

InFile»cf  1 ; 

InFile>>alpha; 

CAMdoubleVector  x(J+l) ,z(K+l) ,dz(K+l) ,nu(K+l) ; 

pi=4.0*atan2(l .0,1.0) ; 

x_init=-27 . 5*pi ; 

x_f inal=27 . 5*pi ; 

z_init=-20.0; 

z_f inal=20.0; 

t_init=0; 

t_f inal=145; 

dx=(x_final  -  x_init)/J; 
dz_eq= (z_f inal-z_init ) /K ; 

mid=long(ceil(K/2.0)) ; 

x(l)=x_init ; 

for  (j=2;  j<=J+l;  j++) 
x(j)=x(j-l)+dx; 

z(l)=z_init ; 

for  (k=2 ;  k<=K+l;  k++) 
z (k) =z (k-1) +dz_eq; 

for  (k=l ;  k<=  K+l;  k++) 

z(k)=20*sinh(z(k) *3/20)/sinh(3) ; 

for  (k=l ;  k<=K;  k++) 
dz(k)=z(k+l)-z(k) ; 


double  fu,  fv,  fw,  gu,  gv,  gw,  u,  v,  w; 
double  vx_init; 
double  II,  12,  13  ,14; 

t ;  ///////////////////////////////////// 

//  Initial  Conditions 
///////////////////////////////////// 
for  (j=l;  j<=J+l;  j++){ 
for  (k=l ;  k<=K+l ;  k++){ 
vx_init=tanh(z(k) ) ; 
v_init (j ,k, 1)=1 .0; 
v_init ( j , k , 2) =vx_init ; 
v_init(j ,k,3)=0; 
v_init(j ,k,4)=0; 
v_init(j ,k,5)=0; 
v_init (j ,k,6)=l .0; 
v_init(j  ,k,7)=0;}}- 

j l=long(ceil (25*pi/dx)+l) ; 
j2=long(ceil(30*pi/dx)+l) ; 

for  (j=j 1 ;  j<=j2;  j++){ 
for (k=l ;  k<=K+l ;  k++) 
v_init(j ,k,2)=v_init (j ,k,2) 

- .008* (sin( . 4*x(j )))*(1.0/(1. 0+pow(z(k) ,2.0))) ;> 

for  (j=l;  j<=J+l;  j++){ 
for  (k=l ;  k<=K+l ;  k++) 

v_init(j ,k,8)=1.0+.5*pow(v_init(j ,k,2) ,2.0) ;} 

/////////////////////////////////////////////////// 

//  time  loop 

/////////////////////////////////////////////////// 

vn=v_init ; 

double  sum_t=0; 

double  dt_cpu=0; 

double  t_start=realtime () ; 

do{ 

for  (j=l;  j<=J+l;  j++){ 
f or (k=l ;  k<=K+l ;  k++){ 
rho_doub=vn( j ,k, 1) ; 
vx=vn ( j , k , 2) / rho_doub ; 
vy=vn ( j , k , 3) / rho_doub ; 
vz=vn ( j , k , 4) / rho_doub ; 

p_doub=vn(j ,k,8)-.5*((pow(vx,2)  +  pow(vy,2) 
+pow(vz,2))/rho_doub)-.5*(pow(vn(j ,k,5) ,2) 
+pow(vn(j ,k,6) ,2)+pow(vn(j ,k,7) ,2)) ; 
A=gamma*p_doub/ rho_doub ; 

B=(pow(vn(j ,k,5) ,2)+pow(vn(j ,k,6) ,2) 

+pow(vn(j ,k,7) ,2) )/rho_doub; 
cf x=sqrt ( . 5* ( A+B+sqrt (pow ( A+B , 2) 

-4*A*pow (vn ( j , k , 5) , 2) / rho_doub) ) ) ; 
cf z=sqrt ( . 5* (A+B+sqrt (pow (A+B , 2) 

-4*A*pow(vn( j ,k,7) ,2)/rho_doub))) ; 
eigx(j ,k)=f abs (vx)+cfx; 
eigz(j ,k)=fabs(vz)+cfz;}} 


dz(K+l)=dz(l) ; 

///////////////////////////////////////// 

//  Solution  Variables 

///////////////////////////////////////// 

CAMdoubleArray  v_init (J+l ,K+1 ,8) ; 

CAMdoubleArray  vn(J+l ,K+1 ,8) ; 

CAMdoubleArray  f (J+l, K+l, 8) ; 

CAMdoubleArray  g(J+l ,K+1 ,8) ; 

CAMdoubleArray  f _p(J+l ,K+1 ,8) ; 

CAMdoubleArray  g_q(J+l ,K+1 ,8) ; 

CAMdoubleArray  v_star ( J+l ,K+1 ,8) ; 

CAMdoubleArray  v_p ( J+l , K+l , 8) ; 

CAMdoubleArray  v_q(J+l ,K+1 ,8) ; 

CAMdoubleArray  f _star ( J+l ,K+1 ,8) ; 

CAMdoubleArray  g_star ( J+l ,K+1 ,8) ; 

CAMdoubleArray  vnpl ( J+l ,K+1 ,8) ; 

double  vx,  vy,  vz,  cfx,  cfz; 

CAMdoubleMatrix  eigx(J+l ,K+1) ,  eigz(J+l ,K+1) ; 

CAMdoubleVector  vn_vector(8) ,fluxx(8) ,fluxz(8) ,v_star_vec(8) ; 


max_eigx=eigx.max() ; 
max_eigz=eigz .max() ; 

dt=cf 1/ ( (max_eigx/dx)+(max_eigz/dz(mid) ) ) ; 

lambda=dt/dx; 

t=t+dt ; 

f or (k=l ;  k<=K+l ;  k++) 
nu(k)=dt/dz(k) ; 

for  (j=l;  j <= J+l ;  j++){ 
for  (k=l ;  k<=K+l ;  k++){ 
for  (1=1;  1<=8;  1++) 

vn_vector (l)=vn(j ,k,l) ; 
f lux_x(vn_vector ,  f luxx) ; 
f lux_z(vn_vector ,  f luxz) ; 

for  (1=1;  1<=8;  l++){ 
f (j ,k,l)=f luxx(l) ; 
g(j  ,k,l)=f luxz(l)  ;»> 

for  (1=1;  1<=8;  l++){ 
for  (k=l ;  k<=K+l ;  k++){ 
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for  (j=2;  j<=J;  j++){ 

f u=alpha* (f ( j +1 , k , 1) -f ( j , k , 1) ) ; 
fv=.5*(f (j+l,k,l)"f (j-l,k,l)) ; 
fw=alpha* (f (j ,k,l)-f (j-l,k,l)) ; 
u=alpha*(vn(j+l,k,l)-vn(j ,k,l)) ; 
v=.5*(vn(j+l,k,l)-vn(j-l,k,l)) ; 
w=alpha* (vn( j ,k,l)-vn(j-l,k,l)) ; 
f_p(j ,k,l)=minmod3(fu,  fv,  fw) ; 
v_p(j ,k,l)=minmod3(u,  v,  w);} 

f_p(l,k,l)=0; 

f_p(J+l,k,l)=0; 

v_p(l,k,l)=0; 

v_p(J+l,k,l)=0;> 

for  ( j =1 ;  j<=J+l;  j++){ 
for  (k=2 ;  k<=K;  k++){ 

gu=alpha* (g( j ,k+l,l)-g(j ,k,l)) ; 
gv=.5*(g(j ,k+l,l)-g(j ,k-l,l)) ; 
gw=alpha* (g( j ,k,l)-g(j ,k-l,l)) ; 
u=alpha* (vn( j ,k+l,l)-vn(j ,k,l)) ; 
v=.5*(vn(j ,k+l,l)-vn(j ,k-l,l)) ; 
w=alpha* (vn( j ,k,l)-vn(j ,k-l,l)) ; 
g_q(j ,k,l)=minmod3(gu,  gv,  gw); 
v_q(j ,k,l)=minmod3(u,  v,  w);} 

g_q(j ,1,1)=0; 
g_q(j,K+l,l)=0; 
v_q( j ,1,1)=0; 
v_q( j  ,K+1,1)=0;» 

/////////////////////////////////////// 

//predicted  values 

/////////////////////////////////////// 
for  (j=l;  j<=J+l;  j++){ 
for  (k=l ;  k<=K+l ;  k++){ 
for (1=1;  1<=8;  l++){ 

v_star(j ,k,l)=vn(j ,k,l)- . 5*lambda*f _p( j ,k,l) 

- . 5*nu(k) *g_q( j ,k,l) ; 
v_star_vec (1) =v_star ( j , k , 1)  ; } 

flux_x(v_star_vec,  f luxx) ; 
flux_z(v_star_vec,  f luxz) ; 

for  (1=1;  1<=8;  l++){ 

f_star(j ,k,l)=f luxx(l) ; 
g_star(j ,k,l)=f luxz(l) ;}}} 

//////////////////////////////////////// 

//corrected  values 

//////////////////////////////////////// 

if (odd){ 

for  (j=l;  j<=J ;  j++H 
for  (k=l ;  k<=K;  k++){ 
for  (1=1;  1<=8;  l++){ 

Il=.0625*(v_p(j ,k,l)-v_p(j+l,k,l)) 
-.5*lambda*(f_star(j+l,k,l)-f_star(j ,k,l)) ; 

12= . 0625* (v_p ( j ,k+l ,l)-v_p(j+l ,k+l ,1)) 
-.5*lambda*(f _star (j+1 ,k+l ,l)-f _star (j ,k+l ,1) ) ; 
13= . 0625* (v_q( j ,k,l)-v_q(j ,k+l,l)) 
-.5*nu(k)*(g_star(j ,k+l ,l)-g_star(j ,k,l)) ; 

I4=. 0625* (v_q( j+1, k,l)"V_q( j+1, k+1,1)) 

- . 5*nu (k) * (g_star ( j  +  1 , k+1 , 1) -g_star ( j+1 , k , 1) ) ; 
vnpl(j ,k,l)=.25*(vn(j+l ,k,l)+vn(j ,k,l) 

+vn(j+l ,k+l ,l)+vn(j , k+1,1))  +I1+I2+I3+I4;}}} 

for  (k=l ;  k<=K;  k++){ 
for (1=1;  1<=8;  l++){ 

13= . 125* (v_q( J+1, k,l)-v_q( J+1, k+1,1)) 
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-nu(k) * (g_star (J+1 ,k+l ,l)-g_star ( J+1 ,k, 1) ) ; 

vnpl (J+1 ,k,l)=. 5*(vn(J+l ,k,l)+vn(J+l ,k+l ,1) )  +13 ;» 

for  ( j =1 ;  j<=J;  j++){ 
for  (1=1;  1<=8;  l++){ 

Il=. 125*(v_p(j , K+1, l)-v_p( j+1, K+1,1)) 
-lambda*(f_star(j+l,K+l,l)-f_star(j , K+1,1)) ; 
vnpl (j ,K+1 ,1)=. 5*(vn(j+l ,K+1 ,l)+vn(j , K+1,1))  +11 ; » 

for  (1=1;  1<=8;  1++) 

vnpl(J+l, K+1, l)=vn(J+l, K+1,1) ; 

vn=vnpl ; 
odd=f alse ; } 

else-C 

for  ( j  =2 ;  j <= J+1 ;  j++){ 
for  (k=2;  k<=K+l;  k++){ 
for  (1=1;  1<=8;  l++){ 

11= . 0625* (v_p ( j -1 , k-1 , 1) -v_p ( j ,k-l,l)) 

-.5*lambda*(f _star (j ,k-l ,l)-f _star (j-1 ,k-l ,1) ) ; 
I2=.0625*(v_p(j-1 ,k,l)-v_p(j ,k,l)) 
-.5*lambda*(f_star(j ,k,l)-f_star(j-l,k,l)) ; 

I3=. 0625* (v_q( j-1, k-1, l)-v_q( j-1, k,l)) 

- .5*nu(k) * (g_star (j-1 ,k, 1) -g_star (j-1 ,k-l ,1) ) ; 

14= . 0625* (v_q( j ,k-l ,l)-v_q(j ,k,l)) 

- . 5*nu(k) * (g_star (j ,k, l)-g_star (j ,k-l,l)) ; 
vnpl(j ,k,l)=.25*(vn(j ,k-l ,l)+vn(j-l ,k-l ,1) 

+vn(j ,k,l)+  vn(j-l ,k,l) )  +  I1+I2+I3+I4; »> 

for  (k=2;  k<=K+l;  k++){ 
for  (1=1;  1<=8;  l++){ 

I3=.125*(v_q(l,k-l,l)-v_q(l,k,l)) 

-nu(k) * (g_star (1 ,k,l)-g_star (1 ,k-l ,1) ) ; 
vnpl  (1  ,k,l)  =  .  5*  (vn(l  ,k-l  ,l)+vn(l  ,k,  1)  )  +13 ; }-} 

for  ( j  =2 ;  j <= J+1 ;  j++){ 
for  (1=1;  1<=8;  l++){ 

11= . 125* (v_p(j-l , 1 ,l)-v_p(j ,1,1)) 

-lambda* (f_star(j , 1 ,l)-f _star (j-1 , 1,1)) ; 
vnpl(j ,l,l)=.5*(vn(j-l,l,l)+vn(j , 1,1)) +I1;}} 

for  (1=1;  1<=8;  1++) 

vnpl (1 , 1 ,l)=vn(l , 1,1); 

vn=vnpl ; 
odd=true ; } 

e_int (s) =0 ; 

for(j=l;  j<=J+l;  j++)f 
f or (k=l ;  k<=K+l ;  k++){ 

e_trans=.5*(pow(vn(j ,k,4) ,2.0)/vn(j ,k, 1) )*dx*dz(k) ; 
e_int (s)=e_int (s)+e_trans ; }} 

time(s)=t ; 
s++; 

dt_cpu=realtime()-t_start ; 
sum_t=sum_t+dt_cpu ; 
t_start=realtime () ; 

}while(t<t_f inal) ; 

writeout (vn,x,z, J,K) ; 

return  0; 

> 

/////////////////////////////////// 
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